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Abstract 

In a progress toward searching for the QCD critical point, we study the finite density phase 
transition of Nj = 4 and 2 lattice QCD at finite temperature with the canonical ensemble approach. 
We develop a winding number expansion method to accurately project out the particle number 
from the fermion determinant which greatly extends the applicable range of baryon number sectors 
to make the study feasible. Our lattice simulation was carried out with the clover fermions and 
improved gauge action. For a given temperature, we calculate the baryon chemical potential from 
the canonical approach to look for the mixed phase as a signal for the first order phase transition. 
In the case of Nf = 4, we observe an "S-shape" structure in the chemical potential-density plane 
due to the surface tension of the mixed phase in a finite volume which is a signal for the first order 
phase transition. We use the Maxwell construction to determine the phase boundaries for three 
temperatures below Tc- The intersecting point of the two extrapolated boundaries turns out to be 
at the expected first order transition point at Tc with fi = 0. This serves as a check for our method 
of identifying the critical point. We also studied the Nj = 2 case, but do not see a signal of the 
mixed phase for temperature as low as 0.83 Tc. 

PACS numbers: ll.15.Ha, ll.30.Rd 
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I. INTRODUCTION 

Quantum Chromodynamics (QCD) is a fundamental theory which describes the strong in- 
teraction of quarks and gluons, from which the nucleons are made of. The knowledge of QCD 
at finite temperature and finite baryon chemical potential is essential to the understanding 
of a variety of phenomena. Exploring the QCD phase diagram, including identification of 
different phases and the determination of the phase transition line is currently one of the 
intensely studied topics [l|, l2(]. Guided by phenomenology and experiments, many candidate 
phase diagrams have been proposed, one of such conjectured phase diagrams [3] is sketched 
in Fig. [U In certain limits, in particular for large temperatures T or large baryon-chemical 
potential fiB, thermodynamics is dominated by short-distance QCD dynamics and the theory 
can be studied analytically. However, the most interesting phenomena lie in regions where 
non-perturbative features of the theory dominate. The only known systematic approach is 
the first-principle calculation via lattice QCD using Monte Carlo simulations. 
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FIG. 1. Conjectured QCD phase diagram 



The thermodynamics of strongly interacting matter has been studied extensively in lattice 
calculations at vanishing baryon chemical potential. Current lattice calculations strongly 
suggest that the transition from the low temperature hadronic phase to the high temperature 
phase is a continuous but rapid crossover happening in a narrow temperature interval around 



T, 



170 Me V 



10| . On the other part of the phase diagram — large baryon chemical 



potential but very low temperature, a number of different model approaches lll-ll8| suggest 
that the transition in this region is strongly first order, although this argument is less 
robust. This first order phase transition is expected to become less pronounced as we lower 



the chemical potential and it should terminate at a second order phase transition point - 
the critical point. 

The search for the QCD critical point has attracted considerable theoretical and exper- 
imental attention recently. The possible existence of this point was introduced sometime 
ago 13|, uJ^. It is apparent that the location of the critical point is a key to understanding 
QCD phase diagram. For experimental search of the critical point, it has been proposed to 
use heavy ion collisions at RHIC 19|, l20|] . The appearance of this point is closely related to 
hadronic fluctuations which may be examined by an event-by-event analysis of experimental 
results. The upcoming RHIC lower energy scan and FAIR (GSI) will focus on the region 
T > 100 MeV and fiB ~ 600 MeV where the critical point is predicted to exist in theoreti- 
cal models. However, in order to extract unambiguous signals for the QCD critical point, 
quantitative calculations from first-principle lattice QCD are indispensable. 

Remarkable progress in lattice simulations at zero baryon chemical potential has been 
made in recent years; however, simulations at non-zero chemical potential are difficult due to 
the complex nature of the fermionic determinant at non-zero chemical potential. The phase 
fluctuations produce the notorious "sign problem" . The majority of current simulations are 
focusing on the small chemical potential region fiq/T <^ 1 where the "sign problem" appears 
to be controllable. Most of them are based on the grand canonical ensemble (T, fiB as 
parameters). Since the existence and location of the critical point is still unknown, we need 
an algorithm which can be applied beyond the small chemical potential region . This is one 
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of the motivations for our studies via the canonical ensemble approach 

We have proposed an algorithm based on the canonical partition function to alleviate the 
overlap and fluctuation problems |26|, l27| . The method we use is computationally expensive 
since every update involves the evaluation of the fermionic determinant; however, finite 
baryon density simulation based on this method has been shown to be feasible at temperature 
above ~ 80% Tc where the sign problem is under control 27|, l28|. In this approach, we 
measure the baryon chemical potential to detect the phase transition. With the aid of the 
winding number expansion technique 29l-l32| . a program was outlined to scan the QCD 
phase diagram in an effort to look for the critical point 33|, |34 |. 

In the present work, we shall consider the cases with four- and two-flavors of degenerate 
quarks. It is numerically easier to carry out HMC simulations with even number of flavors. 
Furthermore, the four-flavor case is known to have a first order phase transition at finite 



temperature with /x = and the phase transition hne is extended to the finite /i. We shall 
use the Maxwell construction to determine the phase boundaries of the mixed phase in the 
canonical ensemble approach at several temperatures lower than T^ and extrapolate the 
boundaries to see if they meet at T^ and /i = 0. This provides us a desirable testing ground 
for our method before we launch a more numerically intensive simulation to search for the 
critical point in the Nf = 3 case. Besides, we can apply the method to the two-fiavor case 
which could be similar to the three-fiavor case and, in this case, will likely have a critical 
point which we can attempt to locate. 

In this paper, we present results based on simulations on 6^ x 4 lattices with clover fermion 
action for two and four fiavors QCD with quark masses which correspond to pion masses at 
~ 700 — 800 MeV. We plot the chemical potential as a function of baryon density. In finite 
volume, due to the non-zero contribution from the surface tension, the first order phase 
transition will be signaled by an "S-shape" structure in this plot. The phase boundaries 
of the coexistent phase can be determined by "Maxwell construction" . We clearly observe 
an S-shape structure in the simulation of the Nf = 4 case, indicating a first order phase 
transition |25| . In our present study, we do not see such a structure in the Nf = 2 case down 
to 0.83 Tc. 

The paper is organized as follows: canonical partition function method is outlined in 
Sec. II, we also present winding number expansion methods (WNEM) and simulation pa- 
rameters in this section. In Sec. Ill, we carry out the measurement of the baryon chemical 
potential as an attempt to determine the phase boundaries in order to scan the phase dia- 
grams with Nf = 4 and A'^^ = 2. We then conclude our study in Sec IV and give an outline 
for the future simulation oi Nf = 3 case. 

II. ALGORITHM 

A. Canonical partition function 

The canonical partition function in lattice QCD can be derived from the fugacity expan- 
sion of the grand canonical partition function, 

Z{V,T,i^) = J2MV,T,k)e'^'/^, (1) 



where k is the net number of quarks (number of quarks minus the number of anti-quarks) and 
Zc is the canonical partition function. We note here that on a finite lattice, the maximum net 
number of quarks is limited by the Pauli exclusion principle. Using the fugacity expansion, 
it can be shown that the canonical partition function can be written as a Fourier transform 
of the grand canonical partition function, 

ZciV,T,k) = ^[ d0e-^^'^Z(\/,T,/i)|^=,<^r, (2) 

after introducing an imaginary chemical potential /i = i(f)T. 

For our study, we will focus on Iwasaki improved gauge action and clover fermion action 



from the following papers |35H37l| . They are defined as 

Z(/3,/€,/x)= fvV{detM{U)ffe-^''^^\ (3) 

M^,y{U) = 6^,y - 6^,yCs^K Y^ a^^Ff,^ ~'^J2[('^~ li)Ui{x) 6^^-^^ + (1 + -fi)W{x - i) 6^_--y 
-K [e-^(l - 74)t/4(x) <5,+a,, + e^il + i,)U,\x - 4) 5,_Q , (5) 



where W^i^^(x) and W^^;^^(a;) are 1 x 1 and 1x2 Wilson loops, F^^ = {ffj_^ - fl^)/{8i), 
where f^i, is the standard clover-shaped combination of gauge links, /3 = G/g"^, Ci = —0.331, 



Co = 1 — 8ci. We adopt nonperturbative 0{a) improved Cs^[38|, Csw = 1 + 0.113(|) + 



0.0158(|)2 + 0.0088(|)3 for Nf = 2 and c,^ = 1 + 0.113(|) + 0.0209(|)2 + 0.0047(|)3 for 
Nf = 4 case, fi = fiqa is the quark chemical potential which is introduced at every temporal 



link. One can perform a change of variables 21| . 



ip{x,X4) -)■ ■?/''(x, X4) = e ^'^''■?/'(f, X4) 

i/j^x, X4) -)■ i/j'{x, X4,) = e^''^il){x, Xa) (6) 

to absorb ^ in the fermionic field. This leaves the chemical potential to reside only on the the 
time-direction links on the last time slice. Now the part of the fermion matrix M(x, y, ^i/T) 
which involves the chemical potential becomes 

M(x, y, /i) = -K [e-'^^'(l - 74)^74 (f, Nt) S^,,nA,,i + ^'''^'(l + I^Wa^x, Nt) S,,,iSy,,N,] hy- 

(7) 



In terms of the imaginary chemical potential, it corresponds to a U{1) phase attached to 
the last time slice 

{U^)u{x) = < (8) 

I Uiy{x) otherwise, 

where fiNt in Eq. ([7]) is replaced with i(j). 

After integrating out the fermionic part in Eq. ([2]), we get an expression 

ZciV,T,k)= fv\Je-^^^^MetkM^f'{U), (9) 

where 

detfcM^/ {U) = — f d0 e-'^^ det M{m, (J); Uff , (10) 

27rio 

is the projected determinant with the fixed net quark number k. 

We shall summarize some of the properties of the canonical ensemble: 

• Since the fermion Dirac matrix with the f/(l) phase on one time slice still satisfies 75 
hermiticity, 

75M([/,0)75 = Mt(f/,0), (11) 

the determinant detM(f/, <^) = detM(f/, 0)* is real. 

• From the charge conjugation property of the Wilson-clover fermion matrix 

M{U, (f)) = C-^M{U\ -<PYC, (12) 

where U^ = U* is the link variable under charge conjugation and C is the gamma 
matrix which has the transformation property C^^'-f^C = —'jj^ and the fact that U 
and W^ configurations has equal weight in the path integral, one can show that 

Zc{V,T,k) = Zc{V,T,-k). (13) 

In other words, the canonical partition function for k quarks is equal to that of k 
anti-quarks. Since the canonical partition function is even in k, Eq. (jH]) becomes 

ZciV,T,k)= f V\Je-^'^^'^RedetkM^f{U), (14) 

and Zciy, T, k) is real, but not necessarily positive. Due to the Fourier transform, 
there can be a sign problem at large quark number and low temperature. Eq. ( fT3|) can 



be similarly derived from considering time reversal transformation. Alternatively, one 
can take the cliarge-75 hermiticity transform (CH) property of the fermion matrix 

M([/, 0) = C-S5M*([/^ -0)75C, (15) 

and the fact that the gauge action in invariant under charge conjugation, i.e. Sg{U) = 
Sgllf^) and that U and If^ configurations has equal weight in the path integral to show 
that the partition function Zc{V,T,k) is real. 

Fermion operators will involve projections in different quark sectors. For example, a 
quark bilinear ipTip = ^ /-^ tpj^Tipf has the following expression for its expectation 
value in the canonical ensemble with k quarks [27| 

X I V^V^e-^f^^^'^'^^Tijj (16) 

^ detfcM^/ "^ ^ '' ^ 

k' 

where we have considered the Fourier transform from Eq. (ITOl) and we define 

TTkTM-^ = — [ d^e-'''^TrTM{U^)-\ (17) 

27ryo 

The summation above runs over all integer values of k but the Fourier coefficients are 
non-zero only for —6Vs < k < 6Kj, with Vs the spatial volume in lattice units. As it 
is clear from the convolution in quark number, sectors with different k will contribute 
to fermionic observables. 

The quark number can be evaluated with the conserved vector charge, 

Nf 

-V'/(x)(l-74)?74(a;)^/(x + t)]. (18) 

{■h)k = ^^ /pUe-^«(^)^ / "d0e-*'^'^detM([/,0)^/(-Ar^)Tr[O], (19) 
Zc{k) J 2,TX Jq 

where O = -k[{1 + 74)f/]M-i(f/, 0) _ (j^ _ ^^)u^M-'^{U, 0)]. From 75 hermiticity, one 
can show that 

Tr[0] = -Tr[0*], (20) 

8 



which means that it is imaginary. On the other hand, from CH transformation, we 
find that {Ji)k is reaL Therefore, 

r-27r 



'\JA)h 



Zcik) 



pUe-Ss(c/). 



2tt 



d(f) i-i) sm{k(j)) det M{U, (py^f {-Nf) Tr[0] 



(21) 



which is odd in k as it should for a charge-odd operator. 
Using the equation above, it is straightforward to show that 



{n)^ 



{JA)k 
1 

Zc{k) 



1 



X5Ue-^s(^)^ / d0e" 



27r 



d 



(det M(f/, 



\Nt 



(22) 



For quark numbers which are multiple of 3, the canonical partition function is invariant 
under the transformation where all the temporal links in a time-slice are multiplied 
with a Z^ element. In this simulation, we shall consider the triality zero sector with 
integral baryon numbers k = Sub by implementing a Z^ hopping in the Hybrid Monte 
Carlo simulation 27|, |39| to avoid being frozen in one Z3 sector. 



To simulate Eq. dSj) dynamically, we can rewrite canonical partition function as 

ZciV,T,k)= /"pUe-^«(^MetfcM^/(f/) 

V\J e-'^'^^MetM^f{U)W{U)a{U), (23) 



where 



and 



W{U) 



a{U) 



RedetkM^f{U)\ 
detM^f{U) 

detkM^f{U) 



(24) 



(25) 



RedetfcM^/([/)| 

Our strategy to generate an ensemble is to employ Metropolis accept/reject method based 
on weight W{U) and fold the phase factor a{U) into the measurements. In short, during the 
simulation, the candidate configuration is "proposed" by the standard Hybrid Monte Carol 
algorithm [40|, |41| and then an accept/reject step is used to correct the probability. 

We note that the accept/reject step is designed to be based on the determinant ratio 



which has been shown to alleviate the fluctuation problem 
acceptance rate. 
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42l | and enhance the 



B. Winding number expansion 

The majority of the time in the simulation is spent on the accept /reject step, specifically 
on computing the determinant of the fermion matrix. On the 6^ x 4 lattice, the matrix has 
10368 X 10368 entries. Although the matrix is very sparse, exact determinant calculation is 



very demanding even on this small lattice. An alternative is to use a noisy estimator [281. |42|. 
In this study we use an exact evaluation of the determinant. One technical challenge has 
to do with the Fourier transform; our original approach was to use an approximation where 
the continuous integration is replaced with a discrete sum, i.e. 

det.M^/(f/) ^ -5^e-^'='^MetM(t/^J^/, 0, = ^. (26) 



It was shown that the errors introduced by this approximation are small 3J] for small 
quark numbers. However, there are two problems with this approach: (i) the computation 
time increases linearly with the net quark number, (ii) although evaluating the determinant 
N times should theoretically allow us to compute projection number as large as A; = A^/2, nu- 
merically we found that for large quark numbers, the Fourier components become too small 
to be evaluated with enough precision, even using double precision fioating point numbers. 
It has been shown [29] that the results of the projected determinant for k larger than 20 
would differ significantly for different choice of A^, which signals a numerical instability. To 
see this problem, we use A^ = 208 to evaluate the fermion determinant and calculate Fourier 
projection using discrete Fourier transform. One would expect | RedetfcM(?7)| to decrease 
exponentially since it is proportional to the free energy. We see that this is indeed the case 
for discrete Fourier transform at small quark number, as shown in Fig. |2](left panel). How- 
ever, as the quark number gets close to 30, the projected determinant calculated using the 
discrete Fourier transform flattens out when it reaches the double precision limit of 10~^^. 
This is the onset of numerical instability. 

This happens because the Fourier coefficient decreases very fast as the quark number 
increases and it quickly becomes comparable to machine precision. The accumulation of 
round-off errors makes it impossible to evaluate the projected determinant. This numeri- 
cal challenge led us to develop a new method which should be free of the above numerical 



problem for quark numbers relevant to the phase transition region in this study 29|. The 



idea of our new method is to first consider the Fourier transform of log det M[U,(f)) instead 
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FIG. 2. Numerical instability of discrete Fourier transform with 208 points (Left). Comparison 
from winding number expansion method and discrete Fourier transform with A^ = 16 evaluations 
(Right). 

of the direct Fourier transform of det M(f/, </>). Using an approximation based on the first 
few components of log det M{U,4>), we can then analytically compute the projected deter- 
minant. The efficacy of the method can be traced to the fact that the Fourier components 
of log det M{U, (f)) are the number of terms in the expansion which characterizes the number 
of quark loops wrapping around the time boundary. They are exponentially smaller with 
increasing winding numbers. This is why we can approximate the exponent of the deter- 
minant very accurately with a few terms which, in turn, allows us to evaluate the Fourier 
components of the determinant precisely. 

To see how this works, we look at the hopping expansion of log det M(f/, 0). We start by 
writing the determinant in terms of the trace log of the quark matrix 



detM(f/,(/)) = exp(logdetM(f/,(/.)) = exp(TrlogM(f/, (^)). 



(27) 



It is well known that Tr log M corresponds to a sum of quark loops. We classify these 
loops in terms of the number of times they wrap around the lattice in the temporal direction. 
Since the determinant is real, the loop expansion is 



TrlogM(f/,0)=^L(f/,, 



(28) 



loops 



MU) + Ee^"'^iy„(^) + e-'-'i'WliU)], 



where n is the net winding number of the quark loops wrapping around the time direction 
and Wn is the weight associated with all these loops with winding number n. Ao{U) is the 
contribution with zero winding number. 



11 



Eq. ( 128|) can also be re-written as 



Trlog M{U,(j))= Ao{U) + [^e*"'^iy„(f/) + e-'^'^W^iU)] 

n 

= Ao{U) + J2 ^n cos(n0 + (5„), 



(29) 



where An = 2\Wn\ and 5„ = arg(iy„) are independent of (p. The pictorial contribution of 
the first few An{U) are demonstrated in Fig. [31 Using Eq. (127]) and Eq. fl29l) we get 



det M{U, 0) = e 



Ao+Ai cos{(j}+5i)+A2 cos(24>+S2)+. 



(30) 
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FIG. 3. Demonstration of first few winding loops in the winding number expansion. The arrows 
at the top of each plot show where the U(l) phases reside. The top two plots show quark loops 
which do not pick up the U{1) phase, so they only contribute to the cj) independent term Aq. The 
lower left and lower right plots show quark loops wrapping around the lattice once and twice in 
the time direction which will contribute to Ai and A2 respectively. 

The Fourier transform of the first order in the expansion can now be computed analyti- 
cally. 



2lT 



2tt 



-ikij) Ao+Aicos{4>+Si) __ Ao+ikSi T /a \ 



(31) 



where Ik is the Bessel function of the first kind with rank k. 
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For higher orders in the winding number expansion, we compute the Fourier transform 
based on the truncated Taylor expansion 



27r 

2tt ji oo oo .^ 



Coolfc(^l) + C+oiIa:+i(^i) + C_oiIfc-l(Ai) + C+02Ia:+2(^i) + •■■ (32) 



The projected determinant is written in terms of the hnear combination of Bessel functions, 
the coefficients c can be easily computed analytically. Using Eq. fl32|) and the recursion 
relation for the Bessel function, lk-i{A) = ^lk{A) +1^+1 (A), the winding number expansion 
method (WNEM) can be extended to higher orders. 

For the purpose of generating configurations we use an approximation where we keep 
only 6 winding loops in the WNEM expansion for TrlogM. To compute the projected 
determinant we use the above Taylor expansion where we keep the first 10 orders in Taylor 
expansion for the second and third term in WNEM and only the first order for the higher 
winding number terms. Referring to Eq. flS2]) this means that we keep only terms with 
n < 6, 1712^3 < 10 and rui^^^Q < 1. We tuned these parameters by comparing the results of 
this approximation with a high order approximation (which we regard as exact to machine 
precision) on a set of gauge configurations that were generated in a previous study. In 
the right panel of Fig. |2] we pick a particular configuration and compare the results of this 
approximation, labeled as WNEM 16, with the high order approximation, labeled WNEM 
208. It can be seen that in the range of quark number of interest {k < 40) the approximation 
is very accurate. This approximation only requires 16 evaluations of the determinant for 
each accept/reject step; we will denote the value of the projected determinant computed 
this way with det^M. 

This approximation is successful because the higher order terms contribute little to the 
value of the projected determinant. It is possible that the approximation is less accurate on 
configurations generated using a different set of parametes than the ones we used to tune the 
approximation. To address this concern, we correct the possible errors when we compute 
the observables. All we have to do is to use a slightly different reweighting factor; we replace 
the factor in Eq. fl25l) with 

= <i-^'m (33) 

|RedetfcM(?7)| ^ ^ 
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This removes any bias the approximation might introduce. The only possible issue is to 
make sure that the reweighting phase doesn't introduce too much noise in the observable. 
This can be gauged by determining how much the average ratio | RedetfcM|/| RedetfcM| 
differs from 1. In our simulations we find that this ratio always stays close to one. 

To compute the projected determinant det^M to machine precision we use a high order 
WNEM approximation and numerical integration of the Fourier transform. We note that 
the numerical integration requires using precision higher than machine precision; we use 
GMP library and perform the Fourier transform using 512-bit precision. The high order 
WNEM approximation involves computing the determinant for many phases (64-128) for 
each configuration in our ensembles. Ideally, we would like to use this method to generate 
configurations but this would cost 4-8 times more computational resources. On the other 
hand, to correct the approximation bias we only have to evaluate this on the saved configu- 
rations which represent a small fraction of the configurations proposed during the dynamical 
evolution. 

Before we conclude this section, we would like to compare the merits of WNEM with those 
of discrete Fourier transform. We compute the values of the projected determinant using the 
discrete Fourier transform with A^ = 16 in Eq. ( I26l) . In Fig. |2] (right panel), these results are 
compared with those from a WNEM approximation using the same number of determinant 
evaluations. We see that the discrete Fourier transform is only valid up to fc = 6. This is 
to be compared to WNEM (16) which takes the same computational time and yet does not 
suffer from this problem and, as expected, the evaluated projected determinant continues to 
decrease as the quark number is increased. It is this projected determinant evaluated from 
WNEM that allows us to scan a wide range of densities in the QCD phase diagram. 

C. Simulation parameters 

A major part of the time in these simulations is spent on computing the fermion matrix 



determinant. We compute the determinant using a numerical package, SuperLU [43j, which 
speeds up sparse matrix determinant calculation by a factor of 10 when compared to the 
usual dense matrix routine from LAPACK. In this study, we work on 6^ x 4 lattices with 
quark masses which correspond to pion masses at 700 — 800 MeV. Since the theoretical 
complexity of this algorithm is between Oln^) and 0{n^), the exact determinant calculation 

14 



will be very time consuming for larger matrices. Even with the SuperLU package, the high 
computational cost prevents us from investigating bigger lattices at the present stage. 

The Fourier transform involves the evaluation of the fermion matrix determinant N times. 
Thus the computation cost increase linearly with N. With the aid the WNEM, A^ = 16 
is found to provide precise enough results for current settings. For this study we used a 
machine which has 16 cores per node. We optimize our algorithm to distribute the 16 
determinant calculations with one determinant calculation per core, so that the calculation of 
each determinant is done in a parallel fashion. Comparing to multi-core parallel computation 
of 16 determinants sequentially, we save almost 40% of the computer time. 

We fixed the scale by using Tq [44] and the pion mass is determined using dynamically 
generated ensembles at zero temperature on 12'^ x 24 lattices for different /3. To locate 
the pseudo critical temperature Tc at zero chemical potential, we varied /3 to look for the 
peak of the Polyakov loop susceptibility. We fixed k at an intermediate quark mass value 



45- 



(k = 0.1371) in order to avoid the unphysical parity-broken phase, the"Aoki phase 
SOj. Since our volume in lattice units is small and quark mass is heavy, we present our 
results in terms of the temperature ratio T/T^ rather than converting them to physical 
units. The temperature T/Tc ratio is determined by measuring the lattice spacing a(/3, k) 
at zero-temperature as 



y(^,^) 



Nt X a{j3c, k) 
Nt X a(/3, k) 



(34) 



TABLE I. Simulation parameters for Nf = 2 



/3 


o(fm) 


m^(MeV) 


y-i(fm-=^) 


T/T, 


1.75 


0.321(2) 


710(10) 


0.139(6) 


0.83(1) 


1.77 


0.306(1) 


744(8) 


0.162(4) 


0.87(2) 


1.79 


0.291(2) 


744(11) 


0.188(5) 


0.92(1) 


1.81 


0.270(3) 


719(14) 


0.235(19) 


1 



We list the simulation parameters (/3, lattice spacing a, r/ivr, inverse spatial volume, and 
T/Tc) in Table [I for N} = 2 and Table |TT] for Nf = 4. 

From these tables, we note that the pion mass varies very little with /3, consequently the 
quark mass is roughly the same in all runs. We also note that the quark mass is quite heavy, 
around the strange quark mass. We chose simulation temperatures below the pseudo-critical 
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TABLE II. Simulation parameters for Nf = 4 



/? 


a(fm) 


m^(MeV) 


y-i(fm-=^) 


T/T, 


1.56 


0.351(2) 


831(10) 


0.107(7) 


0.89(1) 


1.58 


0.340(5) 


828(10) 


0.118(5) 


0.91(2) 


1.60 


0.328(2) 


834(8) 


0.132(12) 


0.95(2) 


1.62 


0.312(4) 


848(6) 


0.152(14) 


1 



temperature because the phase transition at finite chemical potential is believed to appear 
below the transition point at zero temperature. 

To determine the location of the phase transition at non-zero baryon density, we pick three 
temperatures and scan the baryon density while monitoring the chemical potential. The "S- 
shape" wiggle in the chemical potential plot signals the presence of the phase transition. We 
vary the density by changing the net quark number and scan the density range between 4 
and 20 times the nuclear matter density. 

For the HMC proposal step, we adopted the algorithm 40| made exact by an ac- 



cept/reject step at the end of each trajectory[4l|. For the updating process, we set the 
length of the trajectories to 0.5 with 6t = 0.01. The HMC acceptance rate was very close 
to 1 since the step length is small. We adjust the number of HMC trajectories between two 
consecutive finite density accept/reject steps so that the acceptance rate at different k stays 
in the range of 15% to 30%. 



III. RESULTS 



A. Sign problem 



As a first check, we examine the seriousness of the potential sign problem in the canonical 
approach. The average sign in Eq. fl33l) is plotted in Fig. HI For Nj = 4 (left panel), we see 
that they range mostly from 0.6 to 0.2. Except for the point at T = 0.89 Tc and ub = 9 
which has a 1.5 sigma away from zero, the others all have more than 3 sigmas above zero. 
Again, for Nf = 2 (right panel), we see that except for the case at the higher ns and lower 
temperature, the average sign for most of the cases is more than three sigmas above zero. In 
view of the fact that the results we are going to show agree with those from the imaginary 
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chemical potential study and the extrapolated boundaries at finite density on Nf = 4 meet 
at the expected Tc at /i = 0, we believe that the sign fluctuation is not a problem. 
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FIG. 4. Left: Average sign of A'^^ = 4 in terms of ns- Right: Average sign of Nf = 2 in terms of 
B. Baryon chemical potential 



Before starting the discussion on the QCD phase diagram, we will flrst present the "tool" 
we used to determine the phase transition in the canonical ensemble: the baryon chemical 
potential. We measure the chemical potential and plot it as a function of the net quark 
number k. Due to the contribution of the surface tension to the free energy at finite volume, 
the chemical potential in the mixed phase region will display an "S-shape" structure. We 
will scan the baryon number at fixed temperature to look for this S-shape wiggle in chemical 
potential which is a signal for a first order phase transition. 

In the canonical ensemble, the baryon chemical potential is calculated by taking the 
difference of the free energy after adding one baryon 27|, i.e. 



(/^)„, 



F{nB + 1) - F{nB) 



where 



{ub + 1) -Ub 

a{U) 

7(f/) 



1 ^^ Zci^B + 3) 
/3 ^ Zci^riB) 



1 1^ (7(^))o 
/3 («(f/))o 



Redet3n^M"/([/) 

|Redet3„3M"/(f/)|' 

Redet3ns+3M"/([/) 



and 



(35) 

(36) 

(37) 



|Redet3„^M"/(?7)| ' 

are measured in the ensemble with ub baryon number where ()^ in Eq. ( 135|) stands for the 
average over the ensemble generated with the measure |Redet3„^M"'-^(^)|. 
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FIG. 5. Schematic phase diagram of four and two flavors in canonical ensemble. 



C. QCD phase diagram 



Before we present our results, we would like to point out the difference between the 
phase diagram in the grand canonical ensemble and the one in the canonical ensemble. We 
plot the expected canonical ensemble phase diagram in Fig. O In the T - p diagram, the 
first order phase transition line in the grand canonical T — fi diagram becomes a phase 
coexistence region which has two boundaries that separates it from the pure phases. The 
two boundaries will eventually meet at one point. For Nj = 4, this point should be located 
at /is = and the transition temperature Tc. For the Nf = 2 and Nf = 3 cases, if there are 
first order phase transitions, this point will be the critical point at nonzero baryon chemical 
potential. Our method determines the position of the boundaries of the coexistence region at 
the temperatures we use in our scanning; we use an extrapolation to locate the intersection 
point. 

The reasons we study the two and four flavors systems are the following: first of all, it 
is easier to simulate an even number of flavors with HMC 411] • Secondly, in the four- flavor 
case, the first order phase transition extends all the way to zero density. We can then test 
our method by extrapolating the non-zero density results back to zero density to check if 
the transition point identified this way agrees with that determined from the the peak of 
the Polyakov loop susceptibility. Moreover, we can compare our results with those from 
simulations using imaginary chemical potential or those relying on a fi/T Taylor expansion; 
these simulations are expected to identify the transition point reliably for small baryon 
densities. Thirdly, the two-flavor case is expected to have a phase diagram similar to that 
of real QCD (see Fig. [5]). 
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For the above reason, we will use the four- and two-flavor simulations as a benchmark to 
demonstrate the methodology of determining the phase boundaries and locating the critical 
point before tackling the more realistic case. 

We start by determining the location of the transition to the quark-gluon plasma phase 
at zero chemical potential. We determine the order of the phase transition as well as the 
location of T^ using the Polyakov loop susceptibility: 

XP = VNt{{P-{P)f) (38) 

In a finite volume, the susceptibilities are analytic functions, even in the regime where 
phase transitions occur. In the infinite volume limit, on the other hand, phase transitions 
reveal themselves through the divergences of susceptibilities; for crossovers, susceptibilities 
remain finite. The order of the transition can be determined by finite size scaling of the 
susceptibilities. The susceptibility at peak Xpeak behaves as Xpeak oc V", with a being the 
critical exponent. If a = 0, the transition is just a crossover; if < a < 1, it is a second order 
phase transition and if a = 1, it is a first order phase transition. We run simulations for 
two different volumes 6^ x 4 and 10^ x 4 and found a = 0.45(3) for Nj = 2 and a = 1.01(3) 
for Nf =4. It is clear that, the Nf = 4 case has a first order phase transition at zero 



chemical potential which is consistent with other calculations at larger volumes Slj as well 



as a theoretical prediction [52 1. 



We scan the phase diagram by fixing the temperature below T^ while varying ub- Once 
we enter the coexistence region in a finite volume, the contribution from the surface tension 
causes the appearance of a "double-well" effective free energy which leads to an S-shaped 



behavior in the chemical potential vs. baryon number plot [53|. In the thermodynamic limit, 
the surface tension contribution goes away since it is a surface term and the free energy scales 
with the volume; the chemical potential will then stay constant in the mixed-phase region. 
We sketch the expected Nf = 4 phase diagram and mark the boundaries pi and p2 at a 
temperature below T^ in Fig. O The baryon chemical potential in the thermodynamic limit 
is shown in the insert. 

Our results for the baryon chemical potential for Nf = 4 are presented in Fig. [7] for three 
different temperature below T^. It is clear that the chemical potential exhibits "S-shaped" 
wiggles for ub between 5 and 11. To identify the boundaries of the coexistence region and the 
baryon chemical potential of the coexistence phase, we rely on the "Maxwell construction" : 
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p. p.p 



FIG. 6. Schematic plot illustrating the scanning we use to locate the boundaries of the mixed 
phase for QCD with Nf = 4. The infinite volume expectation for chemical potential as a function 
of density is shown in the insert. 



referring to Fig. [7] and Fig. [HI the coexistence chemical potential JIb is the one that produces 
equal areas between the curve of the chemical potential /i^ as a function of ub and the 
constant /Tb line which intersects with fiB at tib, and n^j. This procedure was first used in 
this context by de Forcrand and Kratochvila J25 |. 

Since we do not have a functional description of how the chemical potential depends on 
the density as we cross the coexistence region, we select several points in the S-shape region 
and fit them with a third-order polynomial. This approximation produces good results since 
the value of the coexistence chemical potential, JTb, and the boundary points ub^ and nB2 are 
fairly insensitive to the fitting function and the fit region; the simple third order polynomial 
fit is sufficient. 

We perform the Maxwell constructions for the three temperatures we studied for 0.89 
Tc, 0.91 Tc and 0.95 Tc and the results are presented in Fig. [71 Having determined ub^ and 
nB2 for three different temperatures, we plot the boundaries of the coexistence region and 
perform an extrapolation in hb to locate the intersection of the two boundaries. Although 
there is no "critical point" for the Nf = 4 case, it is expected that the two coexistence phase 
boundaries should cross at zero chemical potential and T = T^. To determine the crossing 
point, we fit the boundary lines using a simple even polynomial in baryon density 

Up) 



TM 



a{Nf,m,){Vp)' + 0{{Vp)') 



(39) 



to do the extrapolation. The phase boundaries and their extrapolations are plotted in 
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FIG. 7. Nf = 4 phase scan along with the dashed hne indicating the coexistence chemical potential 
JJ^/T computed using the Maxwell construction. 
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FIG. 8. Maxwell constructions for T = 0.89 Tc with the horizontal dashed line indicating the 
constant Jib/T and the vertical lines indicating the mixed phase boundaries at ub^ and n^j- 

Fig. H We find the intersection point at rc(p)/Tc(0) = 1.01(5) and p = 0.05(10) fm~^ which 
is consistent with our expectation of p = and T = Tc. 

We would like to point out that the critical temperature was determined using a different 
set of simulations. We ran simulations at zero baryon density and monitored the Polyakov 
loop susceptibility as we increased the temperature. The critical temperature was deter- 
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mined using the peak of the Polyakov loop susceptibihty. The fact that the intersection 
point is consistent with the critical temperature determined from a set of different simula- 
tions constitutes a verification for our method of locating the crossing point of the mixed 
phase boundaries. 




Pb [fm ] 



FIG. 9. Phase boundaries in the temperature vs. density plot for Nj = 4. 



We compare our results to those from a study based on an analytic continuation from 
the imaginary chemical potential with staggered fermions S^l- The results are presented in 
Fig.ra(left panel). The solid eurves represent the staggered fermion results with error barrd 

and our results are given in solid circles - we see that the agreemenio is quite good in spite 
of the fact that the quark mass used in the staggered fermion study, m^ ^ 300 — 400 MeV, is 
significantly smaller than ours at ~ 800 MeV. This is not too surprising in view of the fact 



that t 
mass 



le phase transition curve for Nf = 4 is known to be fairly independent of the quark 



56| 



We also compare our results to the ones from a canonical ensemble study that uses 
staggered fermions J25| (Fig.[TOl right panel). We find that our coexistence region is narrower. 
The upper boundaries between the mixed phase and the quark-gluon plasma phase seem to 
be consistent, but our lower boundary between the hadron phase and the mixed phase lies 
higher in baryon density than that found in the staggered fermion study. The discrepancy 
could come from the difference in the quark masses. The quark mass in our study corresponds 
to ~ 800 MeV pion mass, while for the staggered case it corresponds to rrin ~ 300 MeV . 



^ We also note the agreement with a recent study of Nf = 4 via imginary chemical potential approach 55 |. 
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FIG. 11. Baryon chemical potential for Nf = 2. 

For the Nf = 2 case, the phase diagram (see Fig. [5]) is expected to be similar to the 
conjectured real QCD phase diagram. For certain nonzero quark mass, it also possesses a 
crossover behavior at vanishing chemical potential, a critical point which ends the first order 
phase transition line. For this system, we scan three temperatures below Tc at 0.83 Tc, 0.87 
Tc, and 0.92 Tc. The results are presented in Fig. [TTl We do not observe a clear signal for 
the S-shape structure within the scanned temperature and density range. This may be due 
to the fact that the transition is milder than the sensitivity of our data, or it may be due to 
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the fact that the transition occurs at lower temperatures than the ones we studied - there 
is at least one study that claims that the critical point for Nf = 2 occurs at temperatures 



below 0.8 Tc |53|. 



IV. CONCLUSION AND OUTLOOK 

In this paper, we demonstrate that the canonical ensemble can be used to investigate QCD 
phase diagram at finite temperature and nonzero baryon chemical potential. The procedure 
we implemented to identify the first order phase transition relies on scanning in baryon 
number at a fixed temperature and volume in order to search for the unique mixed phase 
signal in a finite volume as illustrated by an "S-shaped" structure. The well-established 
Maxwell construction is then utilized to determine the boundaries between the mixed phase 
and the pure phases. Finally, the boundaries so determined from several temperatures are 
extrapolated in baryon number and temperature to locate the intersecting point. 

In order to scan the QCD phase diagram at large baryon number, we have developed 
a winding number expansion method to calculate the projected determinant. The winding 
number expansion method expands the scanning region to baryon numbers as large as 12 
in our study region without triggering numerical instability which has rendered the discrete 
Fourier transform impractical. We also note here that the properties of this projection were 



recently investigated by a different group [31 



32|. 



We presented our results for simulations of four and two flavors QCD. In the Nf = 4 
case, the first order phase transition at nonzero chemical potential appears as an "S-shape" 
structure in the plot of chemical potential vs. baryon number (baryon density). The "S- 
shape" structure is related to a finite volume effect which can be understood in terms of 
a surface tension. Although, to confirm the presence of a phase transition, the infinite 
volume limit needs to be studied, we interpret the appearance of the "S-shape" in the finite 
volume as a strong indication of a first order phase transition in the thermodynamic limit. 
Maxwell construction is used to determine the phase boundaries for the Nf = 4 case at three 
temperatures below Tc. 

The intersecting point from the extrapolation of the two boundaries turns out to fall on the 
expected first order phase transition point at T = T^, and zero baryon density. Furthermore, 
our results for the critical chemical potential agree with those from the imaginary chemical 
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potential study with staggered fermions. These facts show convincingly that our method 
can be used to determine the location of the critical point. 

In the Nf = 2 case, we do not see any signal of first order phase transition for temperatures 
as low as 0.83 T^. 

Our present study is carried out on a small volume (6^ x 4) with lattice spacings a = 
0.27 — 0.32 fm and quark masses around that of the strange quark mass. To get closer 
to the physical point, one needs to have lower quark masses, larger volumes and several 
lattice spacing to extrapolate to the continuum limit. This is very challenging since larger 
volumes will require larger baryon numbers for the same density and this is likely to increase 
the severity of the sign problem and lower the Monte Carlo acceptance rate. On the other 
hand, larger volumes will allow us to use smaller quark masses which, in turn, may lower 
the baryon density for the mixed phase due to larger baryon sizes. Noise estimator for the 
fermion determinant 28| or better algorithms may be needed to meet this challenge. Our 
next goal is to extend our present study to the Nf = 3 case. Even with small volumes 
and large quark masses, it is interesting to check whether there exists a first order phase 
transition for this system and, if it exists, to locate its critical point. 
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